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Abstract 

We develop a computational method based on an Eulerian field called the "reference map" , which 
relates the current location of a material point to its initial. The reference map can be discretized 
to permit finite-difference simulation of large solid- like deformations, in any dimension, of the type 
that might otherwise require the finite element method. The generality of the method stems from its 
ability to easily compute kinematic quantities essential to solid mechanics, that are elusive outside 
of Lagrangian frame. This introductory work focuses on large-strain, hyperelastic materials. After 
a brief review of hyperelasticity, a discretization of the method is presented and some non-trivial 
elastic deformations are simulated, both static and dynamic. The method's accuracy is directly 
verified against known analytical solutions. 



1. Introduction 

Two perspectives are used for modeling solid and amorphous materials in the continuum limit. 
Physicists and rheologists tend to prefer an Eulerian treatment where the material can be viewed 
essentially as a fluid equipped with state variables that provide "memory" (see 0, ) ■ On the other 
hand, many in the engineering community prefer a Lagrangian description, where the material body 
is decomposed as a network of volume elements that deform under stress (see [3, 0]). 

Continuum mechanical laws derive from point-particle mechanics applied to a continuum ele- 
ment. The primitive form is resultantly Lagrangian, though an Eulerian conversion can always be 
asserted — one rewrites the constitutive laws in rate form and expands all material time derivatives 
in terms of fixed-space derivatives. Be that as it may, a rigorous Eulerian switch can be a painstak- 
ing mathematical task. This is especially true of solid-like constitutive laws, which often depend 
on nonlinear tensor operations and coupled history-dependent state variables, leading to unduly 
complicated Eulerian rate expansions [7[ . 

To dodge these difficulties, those preferring Eulerian-frame have generally resorted to approxi- 
mations or added conditions that simplify the final constitutive form. While sometimes warranted, 
the connection back to Lagrangian mechanics becomes clouded, complicating the process of deriving 
physically motivated constitutive behavior. 

In this paper, a field we call the "reference map" is utilized to construct and implement solid- 
like constitutive laws in Eulerian-frame with no added approximations. The way the map provides 
"memory" to the system admits immediate computation of kinematic variables crucial to Lagrangian 
solid mechanics. To maintain a clear presentation, several avenues of motivation are first provided 
that discuss the necessary laws of continuum mechanics and the basic quantities of solid kinematics. 
The theory of large-strain elasticity, hyperelasticity, is then sketched. In particular, by enabling 
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quick access to the deformation gradient tensor, the reference map can be used to accurately compute 
solid deformations without the approximations, ambiguities, or pre-conditions of other Eulerian 
approaches. Three non-trivial deformations are then simulated to verify these points. 



2. Basics 

In Eulerian frame, the flow or deformation of a canonical continuous material can be calculated 
by solving a system of equations that includes: 

Pt + V • (pv) = (1) 
V • T + pg = p(v t + v • Vv) (2) 
T T = T. (3) 

The first equation upholds mass conservation, and the next two, respectively, uphold conservation 
of linear and angular momentum. The flow is described by the velocity field v(x, t) and the stresses 
by the Cauchy stress tensor T(x, t), which includes pressure contributions. A consitutive law is 
then asserted to close the system of equations. 



3. Solid kinematics 

We ultimately intend our approach to apply to any material with a "solid-like" constitutive law. 
By solid-like, we mean specifically laws that express the stress tensor in terms of some kinematic 
quantity that measures the local deformation from some nearest relaxed state. This trait reflects 
the microscopic basis of solid stress as arising from potential energy interactions between material 
microconstituents. The simplest solid-like response is isothermal elasticity, where total deflection 
under loading immediately determines the stresses within [H ■ A less basic example would be elasto- 
plasticity, where internal stresses derive from a small elastic component of the total strain. Here, 
the nearest relaxed state can differ from the original unstressed state and may depend on evolving 
state parameters, temperature, and/or rate [§]. 

To encompass the broad definition above, a continuum description for solid-like materials neces- 
sitates a rigorous way of tracking local relative displacements over some finite time period. Without 
making any "small displacement" approximations, a general and robust continuum framework calls 
for a kinematic field x known as the motion function. Suppose at time t = 0, that a body of material 
is in an unstressed reference configuration. The body then undergoes a deformation process such 
that at time t, an element of material originally at X has been moved to x. The motion is defined 
by x = %(X, f). We say that the body at time t is in a deformed configuration. 

The motion can be used to define the deformation gradient F, which is of crucial importance in 
continuum solid mechanics: 



a *(x,f) _ d Xi (x, t) 

ax or ij ~ dx 3 



F (X,f) = ^^ or F ij = A ^ • ( 4 ) 



Note that we use V for gradients in x only, and always write gradients in X in derivative form as 
above. As per the chain rule, the F tensor describes local deformation in the following sense: If 
(DC represents some oriented, small material filament in the reference body, then the deformation 
process stretches and rotates the filament to c?x = F dX in the deformed body. Also, the evolution 
of F can be connected back to the velocity gradient via 

F = (Vv) F (5) 

where we use ' for material time derivatives. 

Since dctF > for any physical deformation, the deformation gradient admits a polar decom- 
position F = RU where R is a rotation, and U is a symmetric positive definite "stretch tensor" 
obeying U 2 = F T F = C. 
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4. General finite-strain elasticity 

To demonstrate the use and simplicity of the method, this paper shall focus on one broad class 
of materials: large-strain, 3D, purely elastic solids at constant temperature. A thermodynamically 
valid constitutive form for such materials is derivable with only minimal starting assumptions. 
Known as hyperelasticity theory, it has become the preeminent elasticity formulation in terms of 
physicality and robustness. Though other elasticity formulations exist (e.g. hypoelasticity and other 
stress-rate models) the next section will recall how these are in fact specific limitting approximations 
to hyperelasticity theory. A brief review of hyperelasticity is provided below to establish the key 
results and demonstrate the physical basis of the theory (see [5[ for details). An analysis of more 
complex solid-like behaviors (e.g. elasto-plasticity, hardening, thermal elasticity) is left as future 
work. 

In essence, one seeks a noncommittal 3D extension of ID spring mechanics, where total relative 
length change determines the force in a fashion independent of deformation path. To institute this, 
presume that the Helmholtz free-energy per unit (undeformed) volume ip and Cauchy stress T both 
depend only on the local deformation: 

tp = t/j (F) , T = T (F) (6) 

where "is used to designate constitutive dependences on kinematic quantities. We also assume that 
if no deformation has occurred (i.e. F = 1), then T = 0. 

Some helpful physical principles refine these dependences immensely. We enforce frame-indifference 
by restricting the dependences to account for rotations. Suppose a material element is deformed 
by some amount, and then the deformed element is rotated. By the frame-indifference principle, 
the rotation should not affect the free-energy, and should only cause the stress to co-rotate. This 
ultimately restricts Eqs[S]to 

ip = ip{C) and T = R f (C) R T . (7) 

Next, we enforce non-violation of the second law. A continuum-level expression of the isothermal 
second law of thermodynamics can be written as the dissipation law 

pip - T : D < (8) 

where D = (Vv + (Vv) T )/2 is the deformation rate, familiar from fluid mechanics. Following a 
procedure originally developed by Coleman and Noll [3J , one can prove mathematically that Eqs [7] 
uphold Ineq [5] under all imposable deformations only if: 

T = 2(dctF)- 1 fM^F t . (9) 

Likewise, C = must correspond to a local minimum of ip. Eq [9] along with the zero deformation 
hypotheses compose the theory of hyperelasticity. 

The above argument demonstrates how the assertions of frame-indifference and the second-law 
require that the assumed dependences of Eq [5] refine to the form of Eq O Each valid choice of ip 
gives an elasticity law that could represent a continuous elastic solid. The "deductive approach" 
above has become a frequently used tool in materials theory. 

5. Some past Eulerian attempts 

To use hyperelasticity, or deductive solid modeling in general, the ability to calculate F during 
a deformation is crucial. In Lagrangian-frame, each point is "tagged" by its start point X, so F can 
always be computed by differentiating current location against initial. In Eulerian the problem is 
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more subtle, as knowledge of past material locations must somehow be procured. As suggested in 



11] , F can be directly evolved by expanding the material time derivative in Eq[5l giving 

F f +v-VF=(Vv)F. (10) 

Unfortunately this cannot be used to solve the general boundary value problem. The term VF 
can only be computed adjacent to boundaries if F is prescribed as a boundary condition. To assign 
F at a boundary implies that the derivative of motion in the direction orthogonal to the surface 
can be controlled. In the general boundary value problem, this information is outside the realm of 
applicability; stress tractions and displacements/velocity conditions can be applied at boundaries, 
but how these quantities change orthogonal to the surface arises as part of the deformation solution. 
Another approach that also advects the F tensor directly (more factually the tensor F -1 ) is the 



Eulerian Godunov method of Miller and Colella [1J] . The method solves for elastic or elasto-plastic 



solid deformation by treating the equations as a system of conservation laws with a nonconserva- 
tive form for the advection of F _1 . It is a sophisticated, high-order method and has had success 
representing solid dynamics and deformation, but is aimed primarily at unbounded domains where 
implementation of a boundary condition on F is unneeded. 

For pure elasticity, several Eulerian, rate-based approaches have been developed that avoid 
directly referring to F but add in several approximations/assumptions. Begin by presuming isotropy. 
It can be shown that this reduces Eq[9]to 
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ah \ oh oh I dl 2 



(11) 



where h,h, h & re the principal invariants of the left Cauchy-Green tensor B = FF T [B| . 

One way to uphold Eq 1111 involves first defining a strain measure E = E(B) that, among other 
features, must asymptote to E = 0.5(F + F T ) — 1 in the small displacement limit |F — 1| <C 1. To 
linear order in E, the elasticity law can then be written as 

T = 2GE + A(trE)l = <T(E). (12) 

Taking the material time derivative gives T = ^(E). The chain rule on E generally leads to a long 
expression in terms of F and F, which can ultimately be rewritten as some function of T and Vv. 

Eulerian expansion of T introduces a term v • VT. Once again, the same problem as that 
encountered in Eq [10] occurs; to compute VT, the full Cauchy stress tensor T must be assigned 
at the boundary. While certain components of T can be controlled at a boundary — namely the 
traction vector Tribound — the components describing stresses along a plane orthogonal to the 
surface cannot, in general, be prescribed. 

To dodge this difficulty the term v- VT is presumed to be negligible. As a consequence of neglect- 
ing stress convection, one accepts certain errors in representing dynamic phenomena. Ultimately, 
what remains is an Eulerian constitutive relation for the evolution of T 

T t = ?f(A(T,Vv)) (13) 

where the function A derives from the choice of strain-measure. While D is sometimes called the 
"strain-rate" , we note that it is not the time rate of change of a valid strain measure; the axes of 
J D dt do not rotate with the material. However, E w D in the small displacement limit for all 
strain definitions. Assuming small displacement and small rate of volume change, Eq 1131 reduces to 
a simple form known as hypoelasticity 

T° = <*f(D) (14) 

where T° is an "objective stress rate" equal to T t plus extra terms that depend on the choice of 
strain measure. 

Hypoelasticity can be seen as a specific approximation to a physically derived isotropic hyper- 
elasticity law. Be that as it may, Eq [TJ] is oftentimes asserted as a starting principle by assigning 
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T°, sometimes arbitrarily, from a list of commonly used stress rates — e.g. Jaumann rate, Truesdell 
rate, Green-Naghdi rate (see [l3j for a detailed review) — thereby cutting off the connection to hy- 
perelasticity. In fact, there are infinitely many stress rate expressions upholding frame- indifference 
that qualify as objective hypoelastic rates. 

Rate forms for elasticity require the assumptions and approximations listed herein, which limit 
their applicability. The neglect of stress convection can pay heavy consequences when attempting 
to represent waves or other dynamic phenomena. While Eq[l3]is fairly general for isotropic linearly 
elastic materials, the resulting equations usually require tedious calculation that must be redone if 
the stress/strain relation is changed. Even in the small strain limit, hypoelasticity's presumptions 
of linearity and isotropy poorly represent some common materials. For instance, granular matter is 
nonlinear near zero strain (due to lack of tensile support), and crystalline solids are not isotropic. 
Rate elasticity, if used as a first principle, also offers no physical basis to account for thermodynamics, 
making it troublesome for theories of thcrmalizcd or non-equilibrium materials. 

6. The reference map 

To sidestep these issues, we now describe a new Eulerian approach to solid mechanics. The key 
is to utilize a fixed-grid field that admits a direct computation of F. Define a vector field called the 
reference map £(x, t) by the evolution law: 



This advection law implies that £ never changes for a tracer moving with the flow. Combined with 
the initial condition, the vector £(x, t) indicates where the material occupying x at time t originally 
started. 

By the chain rule, a material filament obeys dX = (V£)efoc. Thus, 



Altogether, EqsEJ 13 [El and [TBI along with the kinematic expression for the density p = po(detF) -1 , 
compose an Eulerian system that solves exactly for hyperelastic deformation. 

In essence, what we are suggesting is to obtain solid stress in a fashion similar to fluid stress. For 
fluids, the (shear) stress is given by D, which is computed from the gradient of v. Here, we advect 
the primitive quantity £ and use its gradient to construct F. The stress is then obtained from F by 
the constitutive law. This approach alleviates many of the complications discussed previously that 
arise when attempting to directly advect a tensorial quantity like F or T. 

In particular, unlike the advection law of Eq[lGj the reference map is easily definable on bound- 
aries provided complete velocity/displacement boundary conditions. That is, if a boundary point 
originally at Xf, is prescribed a displacement bringing it to x;, at time t, then £(x(,,<) = X&. We 
also note that £ is an integral quantity of F and thus a smoother function. We expect this property 
to be of benefit numerically compared to methods that directly advect F or T. 

The notion of a map that records initial locations of material has been defined by others in 
various different contexts. To these authors' knowledge, it has never been used for the purposes 
of solving solid deformation as described above. Koopman et al. J| use an "original coordinate" 
function akin to our reference map in defining a pseudo-concentration method for flow fronts. The 
inverse of x at time t, which is indeed equivalent to the £ field, is also discussed in Belytschko 0] 
for use in finite element analysis. 

7. Implementation 

In this section, we describe the discretization of the above system of equations. Our general 
strategy is to first evaluate T, then update v using Eq [2 and finally evolve £ with [15] Time 
derivatives in Eqs[2j and 1151 are discretized as a simple Euler step, 



£t + v ■ V£ = for £(x, t = 0) = x = X. 



(15) 



F = (V£) 



(16) 



<9t£(x, t) = (£(x, t + At)- £(x, t)) I At. 
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(a) r (b) 

Figure 1: (a) Angular displacement of washer after inner wall rotation. Inset displays the scalar displacement field; 
note high radial symmetry though the scheme uses a fixed cartesian grid, (b) Displacement field as a function of 
final location after an elastic disk (red outline) is stretched into a triangle. Both tests were performed on a lOOxf 00 
grid. 

On a two-dimensional grid, with grid spacing h, the velocity v and reference map £ are located at 
corner points while stresses T are located at cell centers, (i + |, j + |). Thus, away from any 

boundary, we can compute by finite difference d x £ at the mid-point of horizontal grid edges, and 
similarly, d y $, on vertical grid edges, 

We obtain V£ at cell centers by averaging, allowing us to compute the deformation gradient tensor F 
using Eq[Tni and thus B. We now can define stresses at cell centers by specifying the hyperelasticity 
law. We compute d x T at the mid-point of vertical grid edges, and similarly, d y T on horizontal grid 

edges, 

As a result, we obtain V • T at cell corners, where v is stored. Finally, Vv in equation [5J can be 
discretized in the same manner as V£. Additionally, since Eq[l5]is an advection equation, we use 
a WENO discretization for V£ to guarantee stability [12| . 

In order to solve the system with irregular boundaries, we introduce a level set function, (j>[x). We 
define 4>{x) > inside the solid, <fi(x) < outside, thus implicitly representing the domain boundary 
as the zero level set of 4>(x). Choosing <p(x) to be a signed distance function, i.e. ||V0(a;)|| = 1, 
allows us to compute the cut-cell length a, with < a < 1. For example, if the boundary cuts 

a horizontal cell edge, i.e. 4>i e ft • 4>ri Q ht < 0, then a. — n — rrt^r r- Here, we must change the 

discrete derivatives in Eg 1171 Assuming that 4>i e ft = 4>(i.j) > 0, 

dx£{i+\j) = (bright ~ l ah 

where bright is a boundary condition for Other derivatives near boundaries are treated in the 
same manner. 

8. Results 

8.1. Static solutions 

In this section we present two large-deformation numerical tests, both in plane-strain for sim- 
plicity. To rigorously test the method, each case models a non-trivial, inhomogeneous deformation. 
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First, we solve our system in a circular washer geometry for which the outer wall is fixed and the 
inner wall is rotated over a large angle. Under the Levinson-Burgess hyperelasticity law (see below), 
this environment has an analytical solution, which we use to verify the consistency/correctness of 
the method. Second, we solve the deformation of a disk being stretched into a triangular shape, also 
utilizing the Levinson-Burgess law. This has no analytical solution, and demonstrates the method's 
applicability in cases where the reference and deformed boundary sets differ. We focus here on 
static solutions, obtained by enforcing the boundary conditions and waiting for transients to pass. 
Artificial viscosity was added to the stress law to expedite collapse to the static solution. 

The Levinson-Burgess free-energy function, after application of Eq [TTJ induces the following 
stress law under plane-strain conditions Q 

T = A(/ 3 ) B + f 2 (I u I 3 ) 1 (18) 

where I\ = trB, ^3 = detB are invariants of the B tensor, with fi{I^) = G (3 + I//3) /(A\fTj,) 
and f 2 (Ii,h) = G^/h (k/G- 1/6+ (1 -h)/{Alfj) ~ G/3 - n. In the unstrained state, k and G 
represent the bulk and shear moduli. Throughout, we use k = G =100 kPa. 

Circular washer shear. The analytical static displacement field is A6 = A — B/r 2 and Ar = 
where, A and B are constants that fit the boundary conditions A6*;„ = tt/6 and A9 out = 0. The 
graph in Figure [TJa) shows excellent agreement between our numerical solution (sampled along the 
central horizontal cross-section) and the analytical. We have observed equally high agreement levels 
when the inner wall rotation angle is varied. 

Stretched Disk. The unstressed material shape is a disk that is inscribed perfectly within its final 
equilateral triangular shape. For Xf, on the triangle edge, the boundary condition for the final 
deformed body is £(x;,) = r ° x ^+d(x^) C °" w here x cen is the disk center, r the disk radius, and d(xf,) 
the distance between the disk edge and the triangle as measured along the radial segment containing 
Xfr. Hence, each point on the disk edge is moved outward radially to the triangle. The final, static 
displacement field is shown in figure QJb). 

8.2. Dynamic solutions 

In this section, we display the method's ability to accurately track the motion of a large-strain 
compression wave. Recall from section [5j that dynamic correctness is sacrificed in many stress- 
rate models by neglecting stress convection terms. The reference map on the other hand, enables 
high-accuracy representation of elasto-dynamics as shall now be demonstrated. 

In order to check for the ability of the method to handle dynamic situations, we choose to produce 
an analytical solution for an elastic compresssion wave. For this purpose, consider a material obeying 
the following large-strain elasticity law, 

T = k(V - 1) 

for V = \/B the left stretch tensor. The material body is a rectangular slab constrained in the 
thickness direction (i.e. plane-strain conditions). The unstressed material density is uniform and 
has a value po. 

Under this constitutive law, the following £ and v fields give an exact, analytical solution for a 
rightward moving compression wave passing through the slab: 

x • £ (x, y, t) = x + -Erf(a; - ct) (19) 

x • v (x, y, t) = c (l - —^—^ . (20) 

Due to symmetry, the y and z components of both fields do not change from their initial, unstressed 
values. The constant c = Wk/pq is the wave speed. This solution invokes a large-strain deformation 
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Figure 2: (a) Global convergence of the norm of the error in £ (connected dots). We observe a second a order 
rate of convergence, (b) Travelling x • v wave through one full period t = [0, 10]. Comparison between analytical and 
numerical solutions. 

with compressive strain as high as |x ■ (V — l)x| w 36% at the center of the pulse. Consequently, 
this represents a realistic test of the ability of the present approach to tackle dynamic effects as well 
as large-strain deformation. 

The equations are discretized in space as before, but for the time discretization we embed the 
Euler step described above into a standard second order Runge-Kutta scheme The stability 

restriction of this fully explicit scheme is so that At < amin (Ax, Ay) ^/n/po , for some small 
constant a. From this approach we expect second order global convergence. 

In order to verify the convergence rate, we set up a two-dimensional doubly-periodic domain 
£1 = [—5, 5] x [—5, 5]. We use Eqs fT9l and l20l at t = as the initial condition, with c = k = po = 1. 
The travelling wave solution should come back to it original shape and location at t — 10. For a 
sequence of grids with h = Ax — Ay — {j^, Jj, , § }, we set At using a = j^, and compute the 

Looerror of £, i?£ (h) as, 

El(h)= sup \k-£(x,y,t = 10)-k-£(x,y,t = 0)\, 

{x, y }en 

We report in Figure [2ja) the expected second order global convergence. The convergence rate 
between the two finest grids, h = ^ and h = |, is computed to be 1.97. We have also found nearly 
identical convergence properties for the velocity; the x velocity's convergence rate is found to be 
1.99 between the two finest grids. We conclude the scheme is globally second order accurate for all 
dynamic variables, confirming its ability to capture dynamic solutions under large-strain elasticity. 

Finally, Figure H^b) shows one-dimensional cross sections for x • v at different times. The solid 
line represents the exact solution computed from Eq [20] at t = 0, which, by periodicity of the 
domain, also corresponds to the solution at t = 10. We see that the exact solution and numerical 
solution agree well for h ~ a rather coarse grid. We also plot the numerical solution at t = 2 for 
illustrative purposes. 

9. Conclusion and future work 

This work has demonstrated the validity of the reference map for use in reformulating and 
simulating solid deformation under a completely Eulerian framework. There are still several avenues 
of future investigation. Other material models are to be simulated, most notably elasto-plastic laws 
with and without rate-sensitivity. Also, the reference map has potential to simplify the simulation 
of fluid/solid interactions, due to both phases having a similar Eulerian treatment. Our preliminary 
results on this front are promising, and shall be reported in a future paper. Also, a method to 
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institute traction boundary conditions within this framework, especially the traction-free condition, 
would be important future study. This may ultimately be accomplished with a fluid/solid framework, 
by treating traction-free boundaries as surfaces of interaction with a pressure-free, stationary fluid. 
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